Allgemeine Anmerkungen

Gemischte Lineare Modelle

Beispiel 1

# install.packages("readxl")
# library(readxl)

# LMM1 = read_excel("LMM1.xlsx")
LMM1 = read.csv("https://raw.githubusercontent.com/janikasaretzki/CFH_MVV_2025_26/refs/heads/main/Datasets/LMM1.csv")

# View(LMM1)

names(LMM1)
## [1] "ids"        "department" "experience" "salary"
str(LMM1)
## 'data.frame':    100 obs. of  4 variables:
##  $ ids       : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ department: chr  "sociology" "biology" "english" "informatics" ...
##  $ experience: int  0 1 1 8 1 3 8 2 3 2 ...
##  $ salary    : num  36095 55254 59140 78325 74054 ...
unique(LMM1$department)
## [1] "sociology"   "biology"     "english"     "informatics" "statistics"
# lapply(LMM1[,c("department")], unique)

Einfaches lineares Modell

model1 = lm(salary ~ department, LMM1)
# AV (Kriterium): salary, UV (Prädiktor): department

summary(model1)
## 
## Call:
## lm(formula = salary ~ department, data = LMM1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11922.3  -2985.0    -52.7   3150.7  11138.3 
## 
## Coefficients:
##                       Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)              51843       1127  45.983 < 0.0000000000000002 ***
## departmentenglish        10617       1594   6.659        0.00000000178 ***
## departmentinformatics    26080       1594  16.357 < 0.0000000000000002 ***
## departmentsociology      -3826       1594  -2.399               0.0184 *  
## departmentstatistics     29304       1594  18.379 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5042 on 95 degrees of freedom
## Multiple R-squared:  0.8809, Adjusted R-squared:  0.8759 
## F-statistic: 175.6 on 4 and 95 DF,  p-value: < 0.00000000000000022
qqnorm(resid(model1))

shapiro.test(resid(model1))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(model1)
## W = 0.99367, p-value = 0.925

Unconditional Random Effects Model (UREM)

Nachdem das einfache lineare Modell gezeigt hat, dass die Abteilung (department) einen signifikanten Einfluss auf das Gehalt (salary) hat, möchten wir nun die Unterschiede zwischen den Abteilungen flexibler modellieren.

Das UREM erlaubt es, die Abteilungsunterschiede als zufällige Effekte zu betrachten. Es wird nun angenommen, dass jede Abteilung ihr eigenes durchschnittliches Gehaltsniveau hat.

# install.packages("lme4")
# install.packages("lmerTest")

library(lme4)
library(lmerTest)

model2 = lmer(salary ~ 1 + (1|department), LMM1)
summary(model2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: salary ~ 1 + (1 | department)
##    Data: LMM1
## 
## REML criterion at convergence: 1994.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.38295 -0.58402 -0.00795  0.61588  2.22449 
## 
## Random effects:
##  Groups     Name        Variance  Std.Dev.
##  department (Intercept) 221997992 14900   
##  Residual                25422058  5042   
## Number of obs: 100, groups:  department, 5
## 
## Fixed effects:
##             Estimate Std. Error    df t value Pr(>|t|)    
## (Intercept)    64278       6682     4   9.619 0.000653 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Fokus UREM: Varianzkomponbenten, denn: Hieraus lässt sich die ICC (Intra-Klassen-Korrelation) berechnen (der Anteil der Gesamtvarianz, der durch Unterschiede in der Gruppenzugehörigkeit erklärt wird)

Intraklassen-Korrelation (engl. Intraclass-Correlation; ICC)

ICC_model2 = 221997992 / (221997992 + 25422058)
ICC_model2
## [1] 0.8972514

\(\rightarrow\) 90% der Gesamtvarianz der Gehälter wird durch Unterschiede zwischen den Abteilungen geklärt.

Random Intercept Fixed Slope Model

Einfluss individueller Merkmale (z.B. Berufserfahrung) auf die Gehälter

Annahme: Die mittleren Gehälter können je nach Abteilung weiterhin variieren (Random Intercepts) + die Beziehung zwischen Berufserfahrung und Gehalt ist über alle Abteilungen hinweg konstant (Fixed Slope)

model3 = lmer(salary ~ experience + (1|department), LMM1)
summary(model3)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: salary ~ experience + (1 | department)
##    Data: LMM1
## 
## REML criterion at convergence: 1953.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.93140 -0.80670 -0.07917  0.80156  2.12467 
## 
## Random effects:
##  Groups     Name        Variance  Std.Dev.
##  department (Intercept) 215356636 14675   
##  Residual                18880965  4345   
## Number of obs: 100, groups:  department, 5
## 
## Fixed effects:
##              Estimate Std. Error        df t value     Pr(>|t|)    
## (Intercept) 60469.557   6609.551     4.079   9.149     0.000723 ***
## experience    867.468    148.681    94.013   5.834 0.0000000762 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr)
## experience -0.099
# install.packages("ggplot2")
library(ggplot2)

LMM1$random.intercept.fixed.slope = predict(model3)
names(LMM1)
## [1] "ids"                          "department"                  
## [3] "experience"                   "salary"                      
## [5] "random.intercept.fixed.slope"
ggplot(LMM1, aes(x = experience, y = salary, color = department)) +
  geom_point() +
  geom_line(aes(y = random.intercept.fixed.slope)) +
  labs(x = "Experience (years)", y = "Salary", color = "Department"
  ) +
  theme_minimal() +
  theme(
    panel.grid = element_blank(),
    axis.line = element_line()
  )

Random Intercept Random Slope

Realität: Einfluss der Berufserfahrung variiert von Abteilung zu Abteilung?

model4 = lmer(salary ~ experience + (experience|department), LMM1)
summary(model4)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: salary ~ experience + (experience | department)
##    Data: LMM1
## 
## REML criterion at convergence: 1941.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.76050 -0.84375 -0.01161  0.73272  2.32260 
## 
## Random effects:
##  Groups     Name        Variance  Std.Dev. Corr 
##  department (Intercept) 227233505 15074.3       
##             experience     469881   685.5  -0.26
##  Residual                15518494  3939.4       
## Number of obs: 100, groups:  department, 5
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept) 60149.369   6779.438     4.000   8.872 0.000892 ***
## experience    922.801    335.170     3.987   2.753 0.051390 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr)
## experience -0.274
LMM1$random.intercept.random.slope = predict(model4)
names(LMM1)
## [1] "ids"                           "department"                   
## [3] "experience"                    "salary"                       
## [5] "random.intercept.fixed.slope"  "random.intercept.random.slope"
ggplot(LMM1, aes(x = experience, y = salary, color = department)) +
  geom_point() +
  geom_line(aes(y = random.intercept.random.slope)) +
  labs(x = "Experience (years)", y = "Salary", color = "Department"
  ) +
  theme_minimal() +
  theme(
    panel.grid = element_blank(),
    axis.line = element_line()
  )

Chi-Quadrat-Likelihood Ratio Test

# Model 3: Random Intercept, Fixed Slope
# Model 4: Random Intercept, Random Slope

anova(model3, model4)
## Data: LMM1
## Models:
## model3: salary ~ experience + (1 | department)
## model4: salary ~ experience + (experience | department)
##        npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)   
## model3    4 1992.2 2002.6 -992.10    1984.2                        
## model4    6 1986.5 2002.1 -987.23    1974.5 9.7504  2   0.007634 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Beispiel 2

# install.packages("readxl")
# library(readxl)

LMM2_initial = read.csv("https://raw.githubusercontent.com/janikasaretzki/CFH_MVV_2025_26/refs/heads/main/Datasets/LMM2.csv")
names(LMM2_initial)
##  [1] "subject"            "trainer"            "training"          
##  [4] "performance_time0"  "performance_time1"  "performance_time2" 
##  [7] "performance_time3"  "performance_time4"  "performance_time5" 
## [10] "performance_time6"  "performance_time7"  "performance_time8" 
## [13] "performance_time9"  "performance_time10"
str(LMM2_initial)
## 'data.frame':    80 obs. of  14 variables:
##  $ subject           : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ trainer           : int  5 5 5 5 5 5 5 5 5 5 ...
##  $ training          : chr  "control" "control" "control" "control" ...
##  $ performance_time0 : num  38.9 32.9 35.3 39.3 36.8 ...
##  $ performance_time1 : num  37 36.4 36.9 41.7 34.3 ...
##  $ performance_time2 : num  37.9 31.7 37.1 41.4 35.5 ...
##  $ performance_time3 : num  37.6 34 36.4 41.1 33 ...
##  $ performance_time4 : num  38 33.4 36.5 40.3 36.6 ...
##  $ performance_time5 : num  37.2 32.6 34.9 41.1 38.3 ...
##  $ performance_time6 : num  37.1 33.2 38 40.4 35.3 ...
##  $ performance_time7 : num  37.3 35.5 39.6 39.1 33.5 ...
##  $ performance_time8 : num  34.5 34.8 40 36.2 33.7 ...
##  $ performance_time9 : num  37.8 32.3 37.3 39.3 35.3 ...
##  $ performance_time10: num  34.5 32.8 40.9 38.2 35.9 ...
unique(LMM2_initial$training)
## [1] "control"  "training"
head(LMM2_initial)
##   subject trainer training performance_time0 performance_time1
## 1       1       5  control          38.85797          36.99249
## 2       2       5  control          32.89770          36.38478
## 3       3       5  control          35.32570          36.89093
## 4       4       5  control          39.27142          41.73769
## 5       5       5  control          36.79846          34.28143
## 6       6       5  control          43.62853          42.30369
##   performance_time2 performance_time3 performance_time4 performance_time5
## 1          37.92079          37.62518          38.03515          37.22010
## 2          31.68265          33.99295          33.39586          32.62946
## 3          37.07494          36.38965          36.50572          34.93662
## 4          41.39864          41.10501          40.32918          41.05951
## 5          35.51455          32.98684          36.60746          38.27905
## 6          43.47908          42.67744          40.35285          41.61541
##   performance_time6 performance_time7 performance_time8 performance_time9
## 1          37.10293          37.25276          34.47056          37.76395
## 2          33.21147          35.47573          34.79428          32.33741
## 3          38.00423          39.64092          39.97724          37.32301
## 4          40.37384          39.10527          36.20119          39.33430
## 5          35.30488          33.45806          33.70134          35.29044
## 6          40.39173          40.59770          43.90176          40.65782
##   performance_time10
## 1           34.46748
## 2           32.83205
## 3           40.88486
## 4           38.21026
## 5           35.94133
## 6           39.10271
# WIDE-FORMAT >> LONG-FORMAT

# install.packages("tidyverse") 
# library(tidyverse) 

LMM2 = pivot_longer(
  data = LMM2_initial,
  cols = c(performance_time0, performance_time1, performance_time2, 
           performance_time3, performance_time4, performance_time5, 
           performance_time6, performance_time7, performance_time8, 
           performance_time9, performance_time10),    
  names_to = "Time", # Der Name der ursprünglichen Spalten (z. B. "performance_time0") wird in eine neue Variable namens "Time" geschrieben.
  names_prefix = "performance_time", # Der Präfix "performance_time" wird aus den Spaltennamen entfernt, sodass nur die Zeit (z. B. "0", "1", ...) übrig bleibt.
  values_to = "Performance" # Die Werte in den ausgewählten Spalten werden in eine neue Spalte namens "Performance" geschrieben.
)

LMM2$Time = as.numeric(LMM2$Time)

names(LMM2)
## [1] "subject"     "trainer"     "training"    "Time"        "Performance"
head(LMM2)
## # A tibble: 6 × 5
##   subject trainer training  Time Performance
##     <int>   <int> <chr>    <dbl>       <dbl>
## 1       1       5 control      0        38.9
## 2       1       5 control      1        37.0
## 3       1       5 control      2        37.9
## 4       1       5 control      3        37.6
## 5       1       5 control      4        38.0
## 6       1       5 control      5        37.2
# install.packages("writexl")
# library(writexl)

# writexl::write_xlsx(LMM2, "LMM_2.xlsx") # LMM2 im long-Format

Die Daten haben eine Längsschnittstruktur, da für jede Person (subject) die Outcome-Variable (Performance, d.h. Sprintzeit) zu mehreren Zeitpunkten (Time) erfasst wurde.

\(\rightarrow\) Diese Struktur ermöglicht es uns, die individuelle Entwicklung der Sprintleistung im Zeitverlauf zu analysieren.

Zusätzlich enthalten die Daten zwei Gruppierungsvariablen:

  • Trainer

  • Trainingsstatus

Mögliche Fragestellungen: Wie verändert sich die Sprintzeit im Durchschnitt über die Zeit in Abhängigkeit von der Teilnahme an einem intensiven Intervalltraining?

Mögliche Hypothese: Athleten in der Trainingsgruppe zeigen im Durchschnitt eine stärkere Abnahme der Sprintzeit über die Zeit als Athleten in der Kontrollgruppe.

Einfaches Wachstumskurvenmodell (ohne Prädiktor)

Random Intercept Fixed Slope Model

Annahme: Jeder Athlet eine individuelle Ausgangsleistung (Intercept, d.h. anfängliche Sprintzeit) hat, die um einen globalen Mittelwert streut, während die Wachstumsrate (Steigung, d.h. Verbesserung der Sprintzeit) für alle Athleten gleich bleibt.

\(\rightarrow\) Einfacher Einstieg, da es individuelle Unterschiede im Ausgangspunkt (Intercept) berücksichtigt, ohne die Komplexität zusätzlicher Prädiktoren wie Trainingsstatus oder Trainer oder deren Interaktion einzubeziehen. Es erlaubt eine erste Einschätzung der zeitlichen Entwicklung der Sprintzeiten und der Variabilität zwischen den Athleten.

model5 = lmer(Performance ~ Time + (1|subject), LMM2)
summary(model5)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Performance ~ Time + (1 | subject)
##    Data: LMM2
## 
## REML criterion at convergence: 3673.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1338 -0.6552  0.0057  0.6154  3.1002 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 6.683    2.585   
##  Residual             2.804    1.675   
## Number of obs: 880, groups:  subject, 80
## 
## Fixed effects:
##              Estimate Std. Error        df t value            Pr(>|t|)    
## (Intercept)  37.47091    0.30771  94.10418  121.77 <0.0000000000000002 ***
## Time         -0.27717    0.01785 799.00000  -15.53 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr)
## Time -0.290
LMM2$Random.Intercept.Fixed.Slope = predict(model5, re.form = NA)
# Zusatzargument "re.form = NA" - GLOBALE STEIGUNG, keine individuelle STEIGUNG

ggplot(LMM2, aes(x = Time, y = Performance, group = subject)) +
  geom_point(size = 0.1) +  
  geom_line(size = 0.1) +   
  geom_line(aes(x = Time, y = Random.Intercept.Fixed.Slope), color = "black", size = 0.8) +  
  labs(x = "Time", y = "Performance") +
  theme_minimal() +
  theme(
    panel.grid = element_blank(),
    axis.line = element_line(),
    legend.position = "none"
  )

Random Intercept Random Slope Model

model6 = lmer(Performance ~ Time + (Time|subject), LMM2)
summary(model6)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Performance ~ Time + (Time | subject)
##    Data: LMM2
## 
## REML criterion at convergence: 3559.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.15901 -0.64330 -0.02121  0.59911  2.98766 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  subject  (Intercept) 9.13058  3.0217        
##           Time        0.06458  0.2541   -0.52
##  Residual             2.10207  1.4499        
## Number of obs: 880, groups:  subject, 80
## 
## Fixed effects:
##             Estimate Std. Error       df t value             Pr(>|t|)    
## (Intercept) 37.47091    0.34999 78.99981  107.06 < 0.0000000000000002 ***
## Time        -0.27717    0.03234 78.99935   -8.57    0.000000000000674 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr)
## Time -0.547

Konditionales Wachstumskurvenmodell (mit Prädiktor)

Random Intercept Random Slope Model (2 Gruppen)

  • Zusätzlicher Prädiktor: Beeinflusst die Teilnahme an einem Training die Sprintzeit über die Zeit?

Das Modell berücksichtigt:

  • Random Intercept: Jede Person (subject) hat einen individuellen Ausgangswert (Intercept), der um den globalen Mittelwert streut.

  • Random Slope: Jede Person hat eine individuelle Wachstumsrate (Slope), die um die globale Steigung streut.

  • Training als fester Prädiktor: Das Modell testet, ob die Sprintzeit im Zeitverlauf durch die Teilnahme am Training beeinflusst wird. Die Variable Training hat zwei Ausprägungen: Training vs. Control

model7 = lmer(Performance ~ Time * training + (Time|subject), LMM2)
summary(model7)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Performance ~ Time * training + (Time | subject)
##    Data: LMM2
## 
## REML criterion at convergence: 3555.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.13461 -0.65993  0.00125  0.61249  2.96479 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  subject  (Intercept) 9.21415  3.0355        
##           Time        0.05937  0.2437   -0.53
##  Residual             2.10208  1.4499        
## Number of obs: 880, groups:  subject, 80
## 
## Fixed effects:
##                       Estimate Std. Error       df t value             Pr(>|t|)
## (Intercept)           37.26848    0.49707 78.00008  74.977 < 0.0000000000000002
## Time                  -0.19897    0.04430 78.00246  -4.492            0.0000241
## trainingtraining       0.40486    0.70296 78.00008   0.576               0.5663
## Time:trainingtraining -0.15639    0.06264 78.00246  -2.497               0.0146
##                          
## (Intercept)           ***
## Time                  ***
## trainingtraining         
## Time:trainingtraining *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Time   trnngt
## Time        -0.551              
## tranngtrnng -0.707  0.390       
## Tm:trnngtrn  0.390 -0.707 -0.551
LMM2$Random.Intercept.Random.Slope = predict(model7, re.form = NA)

ggplot(LMM2, aes(x = Time, y = Performance, group = subject)) +
  geom_point(aes(color = training), size = 0.1) +
  geom_line(aes(color = training), size = 0.1) +
  geom_line(aes(x = Time, y = Random.Intercept.Random.Slope), color = "black", size = 1) +
  stat_summary(fun = mean, aes(x = Time, y = Performance, group = training, color = training), geom = "line", linetype = "dashed", size = 1) +
  labs(x = "Time", y = "Performance", color = "Group") +
  theme_minimal() +
  theme(panel.grid = element_blank(), axis.line = element_line(), legend.position = "bottom"
  )

Clusteranalyse

Wiederholung: Die Clusteranalyse ist ein exploratives, strukturentdeckendes Verfahren, mit dem Objekte (z.B. Personen) aufgrund von Ähnlichkeit zu voneinander abgrenzbaren Gruppen (sog. Clustern) zusammenfasst werden.

Aufgabe: Es soll untersucht werden, ob sich in einer Gruppe von Personen unterschiedliche Typen identifizieren lassen. Dazu werden ihre Persönlichkeitsprofile entlang der Giant Three (Extraversion, Neurotizismus, Psychotizismus) betrachtet und mittels Clusteranalyse geprüft, ob sich daraus klar unterscheidbare Gruppen bilden lassen.

# CA = read_sav("CA.sav")
CA = read.csv("https://raw.githubusercontent.com/janikasaretzki/CFH_MVV_2025_26/refs/heads/main/Datasets/CA.csv")
names(CA)
## [1] "SEX" "AGE" "P"   "E"   "N"
CA = subset(CA, select = -c(SEX, AGE))
names(CA)
## [1] "P" "E" "N"
CA$P = scale(CA$P)
CA$E = scale(CA$E)
CA$N = scale(CA$N)

z-Transformation (Standardisierung) - warum? Clusteranalysen beruhen auf abstands- bzw. ähnlichkeitsbasierte Verfahren, bei denen die Skalierung der Variablen eine zentrale Rolle spielt. Durch die Standardisierung werden alle Variablen auf eine gemeinsame Skala mit einem Mittelwert von 0 und einer Standardabweichung von 1 gebracht. So wird sichergestellt, dass jede Variable den gleichen Einfluss auf die Distanzberechnung hat (unabhängig von ursprünglicher Einheit, Streuung und Messbereich)

Agglomeratives Hierarchisches Clustering

Hierarchische Struktur von Clustern (d.h. Fokus auf (Un)Ähnlichkeit zwischen Beobachtungen), Bottom-Up-Prinzip: Jedes Objekt wird anfangs als eigenständiges Clusterblatt betrachtet. In jedem weiteren Schritt werden die beiden ähnlichsten Cluster zu einem größeren Cluster zusammengeführt.

Berechnung der (Un-)Ähnlichkeit - Distanzmatrix

(Un-)Ähnlichkeit zwischen jedem Paar von Objekten im Datensatz

d = dist(CA, method = "euclidean") 
# method-Alternativen: "manhattan", "Mahalanobis"

as.matrix(d)[1:4, 1:4]
##          1         2        3         4
## 1 0.000000 1.7360575 2.065952 1.7797259
## 2 1.736057 0.0000000 3.028084 0.8534846
## 3 2.065952 3.0280835 0.000000 2.8067161
## 4 1.779726 0.8534846 2.806716 0.0000000
# 0: zwei Personen sind identisch auf allen Variablen

Es wird eine Distanzmatrix berechnet, welche die paarweisen Distanzen zwischen allen Personen enthält. Die gewählte Methode ist die euklidische Distanz, welche die geometrische Entfernung zwischen Punkten im Merkmalsraum misst.

Durchführung des Clusteralgorithmus

Grundlegende Frage: Wie modelliere ich die Cluster zueinander (Cluster-Distanz-Maß)?

Typischerweise: Complete-Linkage-Methoden (Maximum-Clustering) oder Ward-Methode

hc = hclust(d, method = "ward.D2")
hc
## 
## Call:
## hclust(d = d, method = "ward.D2")
## 
## Cluster method   : ward.D2 
## Distance         : euclidean 
## Number of objects: 299

Nun: Clusterverfahren durchführen! Der Befehl “hclust” erstellt ein Dendrogramm (Baumdiagramm), das die Hierarchie der Gruppenbildung visualisiert. Die Methode “ward.D2” minimiert die Varianz innerhalb der Cluster (berechnet für jedes Clustering die Fehlerquadratsumme und hält diese so gering wie möglich).

# install.packages("factoextra")
library(factoextra)

fviz_dend(hc, cex = 0.2) + coord_flip()

# cex = 0.2: Schriftgröße wird auf 20% der Standardgröße verkleinert (lesbare Darstellung)

In einem Dendorgramm entspricht jedes Blatt einem Datenpunkt. Objekte, die sich ähnlicher sind, werden zu Ästen kombiniert, bis schließlich alle zu einem einzigen Cluster zusammengeführt werden.

Die cophenetische Distanz beschreibt die Dendrogrammhöhe, bei der zwei Objekte erstmals zu einem gemeinsamen Cluster verschmelzen. Je höher dieser Wert, desto unähnlicher sind die Objekte.

coph = cophenetic(hc)
# coph

Modellpassung

Wir können überprüfen, wie gut das Dendrogramm die zugrunde liegenden Daten widerspiegelt, indem wir es mit der ursprünglichen Distanzmatrix vergleichen, die mit der Funktion dist() berechnet wurde.

\(\rightarrow\) Korrelation zwischen cophenetischer Distanz und ursprünglicher Distanzmatrix \(\rightarrow\) Gute Repräsentation bei Werten ≥ .75

cor(d, coph) # .54
## [1] 0.5399856
# hc_complete = hclust(d, method = "complete")
# coph_complete = cophenetic(hc_complete)
# cor(d, coph_complete) # .57
# 
# d_manhattan = dist(CA, method = "manhattan")
# hc = hclust(manhattan, method = "ward.D2")
# cor(d_manhattan, coph) # .53

\(\rightarrow\) Eingeschränkte Interpretierbarkeit der Clusterlösung aufgrund der nur moderaten Übereinstimmung zwischen Dendrogramm und ursprünglicher Distanzmatrix

Nun: Zuweisung der Fälle zu den Clustern: Anzahl der Cluster k bestimmen

Wahl der optimalen Anzahl von Clustern

Die visuelle Inspektion des Dendrogramms reicht zur Bestimmung der optimalen Clusterzahl i.d.R. nicht aus.

  • k bezeichnet die Anzahl der Cluster, die ein Clustering-Algorithmus erzeugt
  • Die Wahl von k beeinflusst die Ergebnisse sowie die Interpretierbarkeit der Clusterlösung maßgeblich
  • Zu wenige Cluster: Heterogene Gruppen; relevante Unterschiede zwischen Objekten werden überdeckt
  • Zu viele Cluster: Übersegmentierung; die resultierenden Cluster sind instabil und schwer interpretierbar

Elbow-Methode

Häufig verwendete Methode: Elbow-Methode \(\rightarrow\) Berechnung der Gesamt-Within-Cluster-Summe der Quadrate (WSS) für verschiedene Werte von k (Maß für die Homogenität innerhalb der Cluster) und Identifikation des „Knicks“ (Elbow), ab dem eine weitere Erhöhung von k nur noch zu einer geringeren Reduktion der WSS führt \(\rightarrow\) Interpretation: Niedrige WSS-Werte sind günstiger, da sie anzeigen, dass die Datenpunkte innerhalb der Cluster enger beieinander liegen und die Cluster somit homogener sind

elbow_plot = fviz_nbclust(x = CA, FUN = hcut, method = c("wss"))
elbow_plot

# Cluster-Tenzend-Plot, um optimale Anzahl der Cluster zu bestimmen
# „FUN = hcut“ bestimmt die Cluster-Methode. „hcut“ steht für hierarchisches Clustering 
# Alternative: „kmeans“
# „method = c(„wss“)“ verwendet Elbow-Methode, basierend auf der WSS

k_optimal = which.max(diff(diff(elbow_plot$data$y))) + 1
elbow_plot + geom_vline(xintercept = k_optimal, linetype = 2, color = "black")

\(\rightarrow\) „Ellbogen“ repräsentiert den Punkt, an dem der Nutzen abnimmt, wenn k weiter erhöht wird

Silhouetten-Analyse

Gängige Alternative (und Ergänzung) zur Elbow-Methode – sie beantwortet die gleiche Grundfrage (optimale Clusteranzahl k), aber mit einer anderen Logik. Die Silhouetten-Analyse bewertet die Kompaktheit innerhalb der Cluster und Trennung zwischen den Clustern (Silhouettenkoeffizient, zwischen -1 und +1).

\(\rightarrow\) Höherer Durchschnittswert = Bessere Clusterlösung \(\rightarrow\) Erlaubt eine direktere Vergleichbarkeit verschiedener k-Werte

fviz_nbclust(x = CA, FUN = hcut, method = c("silhouette"))

\(\rightarrow\) Anzahl der Cluster (k) = 2 (idealer Wert)

Zuweisung der Fälle zu den Clustern

Dendrogramm nach k (= 2) einfärben:

fviz_dend(hc, cex = 0.2, k = 2, color_labels_by_k = TRUE) + coord_flip()

CA$cluster = cutree(hc, k = 2)
names(CA)
## [1] "P"       "E"       "N"       "cluster"

Clusterbeschreibung und Interpretation

# Deskriptivstatistik der Cluster
aggregate(CA[, c("P", "E", "N")], by = list(Cluster = CA$cluster), FUN = mean)
##   Cluster           P          E          N
## 1       1 -0.06465722  0.3159969 -0.4086034
## 2       2  0.20017166 -0.9782917  1.2649914
aggregate(CA[, c("P", "E", "N")], by = list(Cluster = CA$cluster), FUN = sd)
##   Cluster         P         E         N
## 1       1 1.0333320 0.8833778 0.6945988
## 2       2 0.8650689 0.6453664 0.6961974
# Daten werden nach Clusterzugehörigkeit gruppiert
# Deskriptive Statistiken zeigen die charakteristischen Profile der Cluster in den Variablen P, E und N

# ACHTUNG! STANDARDISIERTE WERTE!

CA_raw = read.csv("https://raw.githubusercontent.com/janikasaretzki/CFH_MVV_2025_26/refs/heads/main/Datasets/CA.csv")

# Nun: Clusterzugehörigkeit aus standardisiertem Datensatz auf Rohdaten übertragen

CA_raw$cluster = CA$cluster
# Voraussetzung: Dieselben Fälle in identischer Reihenfolge!

aggregate(CA_raw[, c("P", "E", "N")], by = list(Cluster = CA$cluster), FUN = mean)
##   Cluster        P        E        N
## 1       1 17.36873 23.36431 10.23009
## 2       2 18.91324 15.80365 24.37443
aggregate(CA_raw[, c("P", "E", "N")], by = list(Cluster = CA$cluster), FUN = sd)
##   Cluster        P        E        N
## 1       1 6.026503 5.160297 5.870383
## 2       2 5.045175 3.769941 5.883893

Cluster 1 ist durch höhere E-Werte und niedrigere N-Werte gekennzeichnet, während sich P auf einem ähnlichen Niveau wie das zweite Cluster befindet.

Cluster 2 zeichnet sich durch höhere N-Werte und niedrigere E-Werte aus, bei insgesamt ähnlichen p-Ausprägungen wie Cluster 1.

ggplot(CA, aes(x = as.factor(cluster), y = P)) +
  geom_boxplot() +
  labs(x = "Cluster",
       y = "Psychotizismus") +
  theme(
    panel.background = element_blank(),
    panel.grid = element_blank(),
    axis.line = element_line(color = "black")
  )

ggplot(CA, aes(x = as.factor(cluster), y = E)) +
  geom_boxplot() +
  labs(x = "Cluster",
       y = "Psychotizismus") +
  theme(
    panel.background = element_blank(),
    panel.grid = element_blank(),
    axis.line = element_line(color = "black")
  )

ggplot(CA, aes(x = as.factor(cluster), y = N)) +
  geom_boxplot() +
  labs(x = "Cluster",
       y = "Psychotizismus") +
  theme(
    panel.background = element_blank(),
    panel.grid = element_blank(),
    axis.line = element_line(color = "black")
  )

t.test(P ~ cluster, data = CA)
## 
##  Welch Two Sample t-test
## 
## data:  P by cluster
## t = -2.1641, df = 143.88, p-value = 0.03211
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
##  -0.50671652 -0.02294124
## sample estimates:
## mean in group 1 mean in group 2 
##     -0.06465722      0.20017166
t.test(E ~ cluster, data = CA)
## 
##  Welch Two Sample t-test
## 
## data:  E by cluster
## t = 13.525, df = 166.06, p-value < 0.00000000000000022
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
##  1.105344 1.483233
## sample estimates:
## mean in group 1 mean in group 2 
##       0.3159969      -0.9782917
t.test(N ~ cluster, data = CA)
## 
##  Welch Two Sample t-test
## 
## data:  N by cluster
## t = -17.867, df = 121.72, p-value < 0.00000000000000022
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
##  -1.859032 -1.488158
## sample estimates:
## mean in group 1 mean in group 2 
##      -0.4086034       1.2649914

Explorative Faktorenanalyse

In der Clusteranalyse haben wir versucht, Struktur in Daten zu entdecken, ohne davor zu wissen, wie viele Gruppen es gibt oder wie sie aussehen (explorativ, datengetrieben, keine a-priori Struktur). Die Frage ist: Können wir dieselbe Idee auch auf einer anderen Ebene anwenden?

Persönlichkeitstheorien gehen i.d.R. davon aus, dass viele einzelne Verhaltensweisen und Selbstbeschreibungen auf wenige grundlegende Dimensionen zurückzuführen sind. Im Folgenden Datensatz liegen neun einzelne Skalen, die unterschiedliche Aspekte von Verhalten und Erleben verfassen, vor. Fragestellung: Lassen sich die Zusammenhänge zwischen den neun Skalen auf eine kleinere Anzahl zugrundeliegender, latenter Faktoren zurückführen - und wenn ja, auf wie viele?

EFA = read.csv("https://raw.githubusercontent.com/janikasaretzki/CFH_MVV_2025_26/refs/heads/main/Datasets/EFA.csv")
names(EFA)
##  [1] "SEX"                     "AGE"                    
##  [3] "Risikobereitschaft"      "Impulsivität"           
##  [5] "Verantwortungslosigkeit" "Aktivität"              
##  [7] "Kontaktfreudigkeit"      "Selbstbewusstsein"      
##  [9] "Minderwertigkeit"        "Schwermut"              
## [11] "Besorgtheit"
EFA = subset(EFA, select = -c(SEX, AGE))
names(EFA)
## [1] "Risikobereitschaft"      "Impulsivität"           
## [3] "Verantwortungslosigkeit" "Aktivität"              
## [5] "Kontaktfreudigkeit"      "Selbstbewusstsein"      
## [7] "Minderwertigkeit"        "Schwermut"              
## [9] "Besorgtheit"

Ziel der EFA: Reduktion einer großen Anzahl von Variablen (z.B. Items eines Fragebogens) auf wenige zugrunde liegende Dimensionen (Faktoren). \(\rightarrow\) Vereinfachung der Datenstruktur, ohne wesentliche Informationen zu verlieren

Logik: Items, die hoch miteinander korrelieren, könnten zum selben, dahinterliegenden Konstrukt gehören (“auf diesem laden”). \(\rightarrow\) Mathematisch basiert die EFA also auf Korrelationen zwischen Indikatoren

cor(EFA)
##                         Risikobereitschaft Impulsivität Verantwortungslosigkeit
## Risikobereitschaft              1.00000000   0.47351664              0.44323461
## Impulsivität                    0.47351664   1.00000000              0.39993072
## Verantwortungslosigkeit         0.44323461   0.39993072              1.00000000
## Aktivität                       0.17112859   0.20249461             -0.32450459
## Kontaktfreudigkeit              0.26041822   0.26048243              0.07312967
## Selbstbewusstsein               0.27754924   0.15869634              0.07172468
## Minderwertigkeit               -0.19405554   0.06794166              0.08302704
## Schwermut                      -0.07104343   0.13201378              0.19287894
## Besorgtheit                    -0.17214201   0.20128015              0.11456815
##                          Aktivität Kontaktfreudigkeit Selbstbewusstsein
## Risikobereitschaft       0.1711286         0.26041822        0.27754924
## Impulsivität             0.2024946         0.26048243        0.15869634
## Verantwortungslosigkeit -0.3245046         0.07312967        0.07172468
## Aktivität                1.0000000         0.45186569        0.33243051
## Kontaktfreudigkeit       0.4518657         1.00000000        0.39587257
## Selbstbewusstsein        0.3324305         0.39587257        1.00000000
## Minderwertigkeit        -0.3622716        -0.40195207       -0.47485579
## Schwermut               -0.3474677        -0.31153103       -0.25263399
## Besorgtheit             -0.1665212        -0.24523910       -0.25264940
##                         Minderwertigkeit   Schwermut Besorgtheit
## Risikobereitschaft           -0.19405554 -0.07104343  -0.1721420
## Impulsivität                  0.06794166  0.13201378   0.2012802
## Verantwortungslosigkeit       0.08302704  0.19287894   0.1145681
## Aktivität                    -0.36227155 -0.34746771  -0.1665212
## Kontaktfreudigkeit           -0.40195207 -0.31153103  -0.2452391
## Selbstbewusstsein            -0.47485579 -0.25263399  -0.2526494
## Minderwertigkeit              1.00000000  0.71410876   0.6765718
## Schwermut                     0.71410876  1.00000000   0.7026063
## Besorgtheit                   0.67657185  0.70260633   1.0000000
round(cor(EFA), 2)
##                         Risikobereitschaft Impulsivität Verantwortungslosigkeit
## Risikobereitschaft                    1.00         0.47                    0.44
## Impulsivität                          0.47         1.00                    0.40
## Verantwortungslosigkeit               0.44         0.40                    1.00
## Aktivität                             0.17         0.20                   -0.32
## Kontaktfreudigkeit                    0.26         0.26                    0.07
## Selbstbewusstsein                     0.28         0.16                    0.07
## Minderwertigkeit                     -0.19         0.07                    0.08
## Schwermut                            -0.07         0.13                    0.19
## Besorgtheit                          -0.17         0.20                    0.11
##                         Aktivität Kontaktfreudigkeit Selbstbewusstsein
## Risikobereitschaft           0.17               0.26              0.28
## Impulsivität                 0.20               0.26              0.16
## Verantwortungslosigkeit     -0.32               0.07              0.07
## Aktivität                    1.00               0.45              0.33
## Kontaktfreudigkeit           0.45               1.00              0.40
## Selbstbewusstsein            0.33               0.40              1.00
## Minderwertigkeit            -0.36              -0.40             -0.47
## Schwermut                   -0.35              -0.31             -0.25
## Besorgtheit                 -0.17              -0.25             -0.25
##                         Minderwertigkeit Schwermut Besorgtheit
## Risikobereitschaft                 -0.19     -0.07       -0.17
## Impulsivität                        0.07      0.13        0.20
## Verantwortungslosigkeit             0.08      0.19        0.11
## Aktivität                          -0.36     -0.35       -0.17
## Kontaktfreudigkeit                 -0.40     -0.31       -0.25
## Selbstbewusstsein                  -0.47     -0.25       -0.25
## Minderwertigkeit                    1.00      0.71        0.68
## Schwermut                           0.71      1.00        0.70
## Besorgtheit                         0.68      0.70        1.00
# install.packages("corrplot")
library(corrplot)

corrplot(cor(EFA), method = "color")

# https://colorbrewer2.org/#type=diverging&scheme=PiYG&n=11
corrplot(cor(EFA), method = "color", col = colorRampPalette(c("#c51b7d", "white", "#4d9221"))(200)) 

# Verlauf mit 200 Abstufungen, feine Übergänge zwischen den drei Ankerfarben

Anwendungsvoraussetzungen

Kaiser-Meyer-Olkin (KMO) Kriterium

# install.packages("psych")
library(psych)

# Korrelationen zwischen den Variablen
## Wert > .60 geeignet für EFA

KMO(EFA) # .71
## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = EFA)
## Overall MSA =  0.71
## MSA for each item = 
##      Risikobereitschaft            Impulsivität Verantwortungslosigkeit 
##                    0.65                    0.66                    0.51 
##               Aktivität      Kontaktfreudigkeit       Selbstbewusstsein 
##                    0.59                    0.83                    0.80 
##        Minderwertigkeit               Schwermut             Besorgtheit 
##                    0.79                    0.76                    0.73
# MSA = Measure of Sampling Adequacy (Overall MSA = KMO)

Bartlett-Test auf Sphärizität

cortest.bartlett(cor(EFA)) # p < .001
## $chisq
## [1] 343.2489
## 
## $p.value
## [1] 0.0000000000000000000000000000000000000000000000000008835381
## 
## $df
## [1] 36
# H0: Korrelationsmatrix ist Einheitsmatrix (Variablen korrelieren perfekt MIT sich aber GAR NICHT MIT ALLEN ANDEREN Variablen)

Weitere Voraussetzungen: Stichprobengröße, lineare Beziehungen zwischen den Variablen, Normalverteilung, Intervallskala, keine Multikollinearität (siehe Folien VO)

Bestimmung der Anzahl der Faktoren

Wie viele latente Faktoren sind notwendig, um die Korrelationen zwischen den Variablen zu erklären?

\(\rightarrow\) Zugang über Eigenwerte: Nur Faktoren mit Eigenwert > 1 werden extrahiert Idee: Die Eigenwerte stammen aus der Eigenwertezerlegung der Korrelationsmatrix. Jeder Eigenwert entspricht der erklärten Varianz eines Faktors, d.h. ein Faktor wird nach diesem Kriterium nur dann behalten, wenn er mehr Varianz erklärt als eine einzelne Variable (dann lohnender Faktor)

Kaiser Kriterium

eigen(cor(EFA))
## eigen() decomposition
## $values
## [1] 3.1929551 2.0686210 1.1907609 0.7097815 0.5906573 0.4437106 0.3270208
## [8] 0.2445874 0.2319054
## 
## $vectors
##              [,1]        [,2]        [,3]        [,4]        [,5]        [,6]
##  [1,]  0.20180955 -0.48942333  0.20344713  0.26857384 -0.41978643 -0.57276092
##  [2,]  0.05295321 -0.56382622 -0.19165245  0.30484118 -0.04053450  0.61663804
##  [3,] -0.05434755 -0.50561927  0.48376815 -0.11450019  0.27195504  0.11761254
##  [4,]  0.33517591 -0.02651339 -0.61858851  0.23814739 -0.20414333 -0.07730031
##  [5,]  0.35670852 -0.20780049 -0.26760594 -0.09271811  0.77058157 -0.32603363
##  [6,]  0.34529673 -0.17439342 -0.09410667 -0.80310917 -0.32503994  0.10879561
##  [7,] -0.48333867 -0.11258204 -0.18216456  0.07951369  0.02199774 -0.18949111
##  [8,] -0.44015227 -0.23338639 -0.18078765 -0.26969617 -0.07327041 -0.32769616
##  [9,] -0.40908585 -0.21324391 -0.39857012 -0.18021110  0.03468271  0.10231488
##              [,7]        [,8]        [,9]
##  [1,] -0.04999693 -0.18684477  0.25402497
##  [2,] -0.39873981 -0.04236226 -0.08344843
##  [3,]  0.54950545  0.27218010 -0.17442216
##  [4,]  0.49054791  0.28578525 -0.27861790
##  [5,] -0.19628183 -0.05439266  0.08825754
##  [6,] -0.12730749  0.20320571  0.14672936
##  [7,] -0.19877833  0.72621308  0.33303696
##  [8,] -0.18582194 -0.16396396 -0.68717485
##  [9,]  0.40882395 -0.45718542  0.45608648
eigen(cor(EFA))$values
## [1] 3.1929551 2.0686210 1.1907609 0.7097815 0.5906573 0.4437106 0.3270208
## [8] 0.2445874 0.2319054

Scree Plot

## Eigenwerte werden nach Größe geplottet, man sucht nach Knick ("elbow")
scree(cor(EFA))

Parallel-Analyse

fa.parallel(EFA, fm = "ml")

## Parallel analysis suggests that the number of factors =  3  and the number of components =  3
# Ist beobachteter Eigenwert größer als das, was man durch den Zufall erwarten muss?

Faktorenextraktion und -rotation

# fa(EFA, nfactors = 3, fm = "ml", rotate = "none")$loadings
fa(EFA, nfactors = 3, fm = "ml", rotate = "promax")$loadings
## 
## Loadings:
##                         ML1    ML3    ML2   
## Risikobereitschaft      -0.101  0.672       
## Impulsivität             0.331  0.617  0.266
## Verantwortungslosigkeit         0.864 -0.494
## Aktivität                      -0.230  1.032
## Kontaktfreudigkeit      -0.189  0.209  0.395
## Selbstbewusstsein       -0.280  0.217  0.233
## Minderwertigkeit         0.825              
## Schwermut                0.817  0.162       
## Besorgtheit              0.913         0.176
## 
##                  ML1   ML3   ML2
## SS loadings    2.426 1.752 1.632
## Proportion Var 0.270 0.195 0.181
## Cumulative Var 0.270 0.464 0.646
# "promax": oblique (vs. orthogonale) Rotationsmethode: Faktoren sind nicht unkorreliert, d.h. nicht unabhängig voneinander
  • Faktor ML1: Psychotizismus (Risikobereitschaft: 0.672, Impulsivität: 0.617, Verantwortungslosigkeit: 0.854
  • Faktor ML2: Extraversion (Aktivität: 1.032, Kontakfreudigkeit: 0.395, Selbstbewusstsein: 0.233)
  • Faktor ML3: Neurotizismus (Minderwertigkeit: 0.825, Schwermut: 0.817, Besorgtheit: 0.913)

Kommunalitäten

… sind ein letzter zentraler Baustein, um die Faktorenlösung zu verstehen. Die Kommunalität einer Variable gibt an, wie viel Varianz dieser Variable durch die extrahierten Faktoren gemeinsam erklärt wird, d.h. für jede Variable wird gefragt_ Wie gut wird diese Variable durch die Faktorenlösung erklärt?

  • Hohe Kommunalität: Die Variable ist gut durch die Faktorenstruktur repräsentiert *Niedrige Kommunalität: Die Variable ist schlecht repräsentiert
fa(EFA, nfactors = 3, fm = "ml", rotate = "promax")$communalities
##      Risikobereitschaft            Impulsivität Verantwortungslosigkeit 
##               0.5042159               0.5318089               0.6699846 
##               Aktivität      Kontaktfreudigkeit       Selbstbewusstsein 
##               0.8576953               0.3831535               0.2998385 
##        Minderwertigkeit               Schwermut             Besorgtheit 
##               0.7473298               0.7075320               0.7029236

Konfirmatorische Faktorenanalyse

CFA = read.csv("https://raw.githubusercontent.com/janikasaretzki/CFH_MVV_2025_26/refs/heads/main/Datasets/CFA.csv")

names(CFA)
## [1] "sex"                       "age"                      
## [3] "Arithmetische_Kompetenz"   "Figural_Induktives_Denken"
## [5] "Wortbedeutung"             "Langzeitgedächtnis"

Hypothese: Die Zusammenhänge zwischen den Subtests “Arithmetische Kompetenz”, Figura-Induktives Denken”, “Wortbedeutung” und “Langzeitgedächtnis” lassen sich durch einen gemeinsamen Generalfaktor der Intelligenz erklären (g-Theorie, Spearman).

Modellspezifikation

# install.packages("lavaan")
library(lavaan)

model = "Generalfaktor =~ Arithmetische_Kompetenz + Figural_Induktives_Denken + Wortbedeutung + Langzeitgedächtnis"


fit = cfa(model, data = CFA)
# lavaan schätzt Faktorladungen (wie stark misst jeder Subtest g?) und Fehlervarianzen

summary(fit) # n.s., die Abweichungen zwischen beobachteter und modellimplizierter Kovarianzmatrix ist nicht größer als zufällig erwartet ("Das Modell passt!" :) )
## lavaan 0.6-20 ended normally after 28 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                         8
## 
##   Number of observations                           196
## 
## Model Test User Model:
##                                                       
##   Test statistic                                 2.356
##   Degrees of freedom                                 2
##   P-value (Chi-square)                           0.308
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                               Estimate  Std.Err  z-value  P(>|z|)
##   Generalfaktor =~                                               
##     Arithmetische_Kompetenz      1.000                           
##     Figural_Induktives_Denken    1.231    0.320    3.848    0.000
##     Wortbedeutung                1.110    0.293    3.790    0.000
##     Langzeitgedächtnis           0.965    0.262    3.679    0.000
## 
## Variances:
##                               Estimate  Std.Err  z-value  P(>|z|)
##    .Arithmetische_Kompetenz      1.029    0.125    8.220    0.000
##    .Figural_Induktives_Denken    0.644    0.111    5.776    0.000
##    .Wortbedeutung                0.830    0.114    7.262    0.000
##    .Langzeitgedächtnis           0.787    0.100    7.833    0.000
##     Generalfaktor                0.259    0.105    2.472    0.013
# ACHTUNG! Der Chi-Quadrat Test ist stichprobenabhängig, deshalb betrachtet man zusätzliche Fit-Indizes

Fit-Indizes

summary(fit, fit.measures = TRUE)
## lavaan 0.6-20 ended normally after 28 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                         8
## 
##   Number of observations                           196
## 
## Model Test User Model:
##                                                       
##   Test statistic                                 2.356
##   Degrees of freedom                                 2
##   P-value (Chi-square)                           0.308
## 
## Model Test Baseline Model:
## 
##   Test statistic                                72.216
##   Degrees of freedom                                 6
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.995
##   Tucker-Lewis Index (TLI)                       0.984
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -1122.279
##   Loglikelihood unrestricted model (H1)      -1121.101
##                                                       
##   Akaike (AIC)                                2260.558
##   Bayesian (BIC)                              2286.783
##   Sample-size adjusted Bayesian (SABIC)       2261.440
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.030
##   90 Percent confidence interval - lower         0.000
##   90 Percent confidence interval - upper         0.148
##   P-value H_0: RMSEA <= 0.050                    0.468
##   P-value H_0: RMSEA >= 0.080                    0.345
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.025
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                               Estimate  Std.Err  z-value  P(>|z|)
##   Generalfaktor =~                                               
##     Arithmetische_Kompetenz      1.000                           
##     Figural_Induktives_Denken    1.231    0.320    3.848    0.000
##     Wortbedeutung                1.110    0.293    3.790    0.000
##     Langzeitgedächtnis           0.965    0.262    3.679    0.000
## 
## Variances:
##                               Estimate  Std.Err  z-value  P(>|z|)
##    .Arithmetische_Kompetenz      1.029    0.125    8.220    0.000
##    .Figural_Induktives_Denken    0.644    0.111    5.776    0.000
##    .Wortbedeutung                0.830    0.114    7.262    0.000
##    .Langzeitgedächtnis           0.787    0.100    7.833    0.000
##     Generalfaktor                0.259    0.105    2.472    0.013

\(\rightarrow\) Wenn die drei Indizes CFI, RMSEA und SRMR im akzeptablen Bereich (siehe z.B. Hu, L. T., & Bentler, P. M. (1999). Cutoff criteria for fit indexes in covariance structure analysis: Conventional criteria versus new alternatives. Structural equation modeling: a multidisciplinary journal, 6(1), 1-55. https://doi.org/10.1080/10705519909540118) liegen: Die Generalfaktor-Hypothese wird empirisch gestützt!

Standardisierte Lösung

summary(fit, fit.measures = TRUE, standardized = TRUE)
## lavaan 0.6-20 ended normally after 28 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                         8
## 
##   Number of observations                           196
## 
## Model Test User Model:
##                                                       
##   Test statistic                                 2.356
##   Degrees of freedom                                 2
##   P-value (Chi-square)                           0.308
## 
## Model Test Baseline Model:
## 
##   Test statistic                                72.216
##   Degrees of freedom                                 6
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.995
##   Tucker-Lewis Index (TLI)                       0.984
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -1122.279
##   Loglikelihood unrestricted model (H1)      -1121.101
##                                                       
##   Akaike (AIC)                                2260.558
##   Bayesian (BIC)                              2286.783
##   Sample-size adjusted Bayesian (SABIC)       2261.440
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.030
##   90 Percent confidence interval - lower         0.000
##   90 Percent confidence interval - upper         0.148
##   P-value H_0: RMSEA <= 0.050                    0.468
##   P-value H_0: RMSEA >= 0.080                    0.345
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.025
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                               Estimate  Std.Err  z-value  P(>|z|)   Std.lv
##   Generalfaktor =~                                                        
##     Arithmetische_Kompetenz      1.000                               0.509
##     Figural_Induktives_Denken    1.231    0.320    3.848    0.000    0.627
##     Wortbedeutung                1.110    0.293    3.790    0.000    0.565
##     Langzeitgedächtnis           0.965    0.262    3.679    0.000    0.491
##   Std.all
##          
##     0.449
##     0.616
##     0.527
##     0.484
## 
## Variances:
##                               Estimate  Std.Err  z-value  P(>|z|)   Std.lv
##    .Arithmetische_Kompetenz      1.029    0.125    8.220    0.000    1.029
##    .Figural_Induktives_Denken    0.644    0.111    5.776    0.000    0.644
##    .Wortbedeutung                0.830    0.114    7.262    0.000    0.830
##    .Langzeitgedächtnis           0.787    0.100    7.833    0.000    0.787
##     Generalfaktor                0.259    0.105    2.472    0.013    1.000
##   Std.all
##     0.799
##     0.621
##     0.722
##     0.765
##     1.000
# install.packages("semPlot")
library(semPlot)

semPaths(fit)

semPaths(fit, "std", edge.label.cex = 1, layout = "tree", whatLabels = "std")

  • Standardisierte Faktorenladungen an den Pfeilen
  • Standardisierte Residualvarianzen (Fehlervarianzen)
  • Varianz des latenten Faktors: Skalierung

Streichelte Linie? Referenzladung, Modell wird fixiert (skaliert), irrelevant für standardisierte Faktorladungen

Messinvarianz

names(CFA)
## [1] "sex"                       "age"                      
## [3] "Arithmetische_Kompetenz"   "Figural_Induktives_Denken"
## [5] "Wortbedeutung"             "Langzeitgedächtnis"
# str(CFA$sex)
CFA$sex = as.factor(CFA$sex)

# Drei Stufen der Messinvarianz (heute) im Fokus: Konfigural, Metrisch, Skalare

## Konfigurale Messinvarianz
fit1 = cfa(model, data = CFA, group = "sex")

## Metrische Messinvarianz
fit2 = cfa(model, data = CFA, group = "sex", group.equal = "loadings")

## Skalare Messinvarianz
fit3 = cfa(model, data = CFA, group = "sex", group.equal = c("intercepts", "loadings"))

# Messinvarianzprüfung (Modellvergleich)
lavTestLRT(fit1, fit2, fit3)
## 
## Chi-Squared Difference Test
## 
##      Df    AIC    BIC   Chisq Chisq diff   RMSEA Df diff Pr(>Chisq)  
## fit1  4 2271.6 2350.3  5.0780                                        
## fit2  7 2267.6 2336.4  7.0245     1.9465 0.00000       3    0.58359  
## fit3 10 2271.4 2330.4 16.8817     9.8571 0.15272       3    0.01982 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Strukturgleichungsmodelle

PoliticalDemocracy Datensatz (Bollen, 1989), lavaan Paket: Der Datensatz enthält Informationen zu politischen und wirtschaftlichen Indikatoren in verschiedenen Ländern. Die zentralen Variablen beziehen sich auf die Demokratie und die industrielle Entwicklung im Zeitraum 1965 bis 1969.

Latente Varablen (Faktoren), die im Modell spezifiziert werden sollen: * Demokratie im Jahr 1960 (dem60) * Demokratie im Jahr 1965 (dem65) * Industrielle Entwicklung im Jahr 1960 (ind60)

Diese werden durch manifeste Variablen (Indikatoren) gemessen: * dem60 wird durch die Indikatoren y1, y2, y3 und y4 repräsentiert * dem65 wird durch die Indikatoren y5, y6, y7 und y8 repräsentiert * ind60 wird durch die Indikatoren x1, x2 und x3 repräsentiert

?PoliticalDemocracy

head(PoliticalDemocracy, 10)
##       y1       y2       y3       y4       y5       y6       y7        y8
## 1   2.50 0.000000 3.333333 0.000000 1.250000 0.000000 3.726360  3.333333
## 2   1.25 0.000000 3.333333 0.000000 6.250000 1.100000 6.666666  0.736999
## 3   7.50 8.800000 9.999998 9.199991 8.750000 8.094061 9.999998  8.211809
## 4   8.90 8.800000 9.999998 9.199991 8.907948 8.127979 9.999998  4.615086
## 5  10.00 3.333333 9.999998 6.666666 7.500000 3.333333 9.999998  6.666666
## 6   7.50 3.333333 6.666666 6.666666 6.250000 1.100000 6.666666  0.368500
## 7   7.50 3.333333 6.666666 6.666666 5.000000 2.233333 8.271257  1.485166
## 8   7.50 2.233333 9.999998 1.496333 6.250000 3.333333 9.999998  6.666666
## 9   2.50 3.333333 3.333333 3.333333 6.250000 3.333333 3.333333  3.333333
## 10 10.00 6.666666 9.999998 8.899991 8.750000 6.666666 9.999998 10.000000
##          x1       x2       x3
## 1  4.442651 3.637586 2.557615
## 2  5.384495 5.062595 3.568079
## 3  5.961005 6.255750 5.224433
## 4  6.285998 7.567863 6.267495
## 5  5.863631 6.818924 4.573679
## 6  5.533389 5.135798 3.892270
## 7  5.308268 5.075174 3.316213
## 8  5.347108 4.852030 4.263183
## 9  5.521461 5.241747 4.115168
## 10 5.828946 5.370638 4.446216
names(PoliticalDemocracy)
##  [1] "y1" "y2" "y3" "y4" "y5" "y6" "y7" "y8" "x1" "x2" "x3"

Hypothese: Die industrielle Entwicklung im Jahr 1960 (ind60) beeinflusst sowohl die Demokratie im Jahre 1960 (dem60) als auch die Demokratie im Jahre 1965 (dem65). Es wird zudem angenommen, dass die Demokratie im Jahre 1960 (dem60) einen Einfluss auf die Demokratie im Jahre 1965 (dem65) hat.

Modellspezifikation

Zwei Hauptkomponenten des SEMs: * Messmodell: Beschreibt die Beziehungen zwischen den latenten Variablen und ihren manifesten Indikatoren * Strukturmodell: Beschreibt die Beziehungen zwischen den latenten Variablen

ACHTUNG: In einem SEM können die Residuen (Fehlerterme) zwischen den Indikatoren miteinander korrelieren. Das bedeutet, dass es neben dem Einfluss der latenten Variablen möglicherweise noch andere gemeinsame Ursachen oder Einflüsse gibt, die zu Ähnlichkeiten in den beobachteten Variablen führen. Solche Korrelationen zwischen den Fehlertermen können spezifiziert werden, um diesen zusätzlichen Zusammenhang im Modell zu berücksichtigen. In unserem Beispiel wird angenommen, dass bestimmte Indikatoren innerhalb des Konstrukts “Demokratie” (sowohl für 1960 als auch für 1965) ähnliche, nicht durch das Modell erklärte Einflüsse teilen, weshalb Residualkorrelationen (d.h. Beziehungen zwischen den Fehlertermen) in das Modell integriert werden müssen.

model = "

# Messmodell
ind60 =~ x1 + x2 + x3
dem60 =~ y1 + y2 + y3 + y4
dem65 = ~ y5 + y6+ y7 + y8

# Strukturmodell, Regression
dem60 ~ ind60
dem65 ~ ind60 + dem60

# Residualkorrelationen
y1 ~~ y5
y2 ~~ y4 + y6
y3 ~~ y7
y4 ~~ y8
y6 ~~ y8
"

Standardisierte Lösung

Das Modell kann nun an die Daten angepasst werden. “Das Modell fitten” beschreibt den Prozess, bei dem die Modellparameter so geschätzt werden, dass die Unterschiede zwischen den theoretisch erwarteten und den tatsächlich beobachteten Werten minimiert werden. Dabei wird in der Regel die ML-Methode verwendet, um die besten Schätzungen für die Parameter zu erhalten.

Um das Modell umfassend zu analysieren, verwenden wir in der Zusammenfassung zusätzlich zu den standardisierten Koeffizienten auch die Fit-Indizes und die R-Quadrat-Werte, um die Modellpassung und Erklärungsleistung zu beurteilen.

fit = sem(model, data = PoliticalDemocracy)
summary(fit, fit.measures = TRUE, standardized = TRUE, rsquare = TRUE)
## lavaan 0.6-20 ended normally after 68 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        31
## 
##   Number of observations                            75
## 
## Model Test User Model:
##                                                       
##   Test statistic                                38.125
##   Degrees of freedom                                35
##   P-value (Chi-square)                           0.329
## 
## Model Test Baseline Model:
## 
##   Test statistic                               730.654
##   Degrees of freedom                                55
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.995
##   Tucker-Lewis Index (TLI)                       0.993
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -1547.791
##   Loglikelihood unrestricted model (H1)      -1528.728
##                                                       
##   Akaike (AIC)                                3157.582
##   Bayesian (BIC)                              3229.424
##   Sample-size adjusted Bayesian (SABIC)       3131.720
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.035
##   90 Percent confidence interval - lower         0.000
##   90 Percent confidence interval - upper         0.092
##   P-value H_0: RMSEA <= 0.050                    0.611
##   P-value H_0: RMSEA >= 0.080                    0.114
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.044
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   ind60 =~                                                              
##     x1                1.000                               0.670    0.920
##     x2                2.180    0.139   15.742    0.000    1.460    0.973
##     x3                1.819    0.152   11.967    0.000    1.218    0.872
##   dem60 =~                                                              
##     y1                1.000                               2.223    0.850
##     y2                1.257    0.182    6.889    0.000    2.794    0.717
##     y3                1.058    0.151    6.987    0.000    2.351    0.722
##     y4                1.265    0.145    8.722    0.000    2.812    0.846
##   dem65 =~                                                              
##     y5                1.000                               2.103    0.808
##     y6                1.186    0.169    7.024    0.000    2.493    0.746
##     y7                1.280    0.160    8.002    0.000    2.691    0.824
##     y8                1.266    0.158    8.007    0.000    2.662    0.828
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   dem60 ~                                                               
##     ind60             1.483    0.399    3.715    0.000    0.447    0.447
##   dem65 ~                                                               
##     ind60             0.572    0.221    2.586    0.010    0.182    0.182
##     dem60             0.837    0.098    8.514    0.000    0.885    0.885
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##  .y1 ~~                                                                 
##    .y5                0.624    0.358    1.741    0.082    0.624    0.296
##  .y2 ~~                                                                 
##    .y4                1.313    0.702    1.871    0.061    1.313    0.273
##    .y6                2.153    0.734    2.934    0.003    2.153    0.356
##  .y3 ~~                                                                 
##    .y7                0.795    0.608    1.308    0.191    0.795    0.191
##  .y4 ~~                                                                 
##    .y8                0.348    0.442    0.787    0.431    0.348    0.109
##  .y6 ~~                                                                 
##    .y8                1.356    0.568    2.386    0.017    1.356    0.338
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .x1                0.082    0.019    4.184    0.000    0.082    0.154
##    .x2                0.120    0.070    1.718    0.086    0.120    0.053
##    .x3                0.467    0.090    5.177    0.000    0.467    0.239
##    .y1                1.891    0.444    4.256    0.000    1.891    0.277
##    .y2                7.373    1.374    5.366    0.000    7.373    0.486
##    .y3                5.067    0.952    5.324    0.000    5.067    0.478
##    .y4                3.148    0.739    4.261    0.000    3.148    0.285
##    .y5                2.351    0.480    4.895    0.000    2.351    0.347
##    .y6                4.954    0.914    5.419    0.000    4.954    0.443
##    .y7                3.431    0.713    4.814    0.000    3.431    0.322
##    .y8                3.254    0.695    4.685    0.000    3.254    0.315
##     ind60             0.448    0.087    5.173    0.000    1.000    1.000
##    .dem60             3.956    0.921    4.295    0.000    0.800    0.800
##    .dem65             0.172    0.215    0.803    0.422    0.039    0.039
## 
## R-Square:
##                    Estimate
##     x1                0.846
##     x2                0.947
##     x3                0.761
##     y1                0.723
##     y2                0.514
##     y3                0.522
##     y4                0.715
##     y5                0.653
##     y6                0.557
##     y7                0.678
##     y8                0.685
##     dem60             0.200
##     dem65             0.961
semPaths(fit)

semPaths(fit, 
         what = "std",
         whatLabels = "std",
         edge.color = "black", # Pfade in schwarz
         edge.label.cex = 0.8, # Größe der Labelschrift
         edge.width = 0.2, # Einheitliche Liniendicke
         layout = "tree", # Alternative: circle, spring
         sizeMan = 6, # Größe der manifesten Variablen
         sizeLat = 8, # Größe der latenten Variablen
         )

# what = "std", # Standardisierte Pfade
# whatLabels = "std", # Labels mit standardisierten Werten

Mediation

MED = read.csv("https://raw.githubusercontent.com/janikasaretzki/CFH_MVV_2025_26/refs/heads/main/Datasets/MED.csv")

names(MED)
## [1] "NSB" "ABK" "BDI"
# NSB = Negative Selbstbewertung
# ABK = Abhängigkeitskognition
# BDI = Beck-Depressions-Inventar (Ausprägung der depressiven Symptomatik)

Annahme: Es wird angenommen, dass der Zusammenhang zwischen Abhängigkeitskognitionen und der Ausprägung der Depressivität durch die negative Selbstbewertung vermittelt wird.

\(\rightarrow\) Hypothese im Sinne einer vollständigen Mediation

# https://lavaan.ugent.be/tutorial/mediation.html

model <- ' # direct effect
             BDI ~ c*ABK
           # mediator
             NSB ~ a*ABK
             BDI ~ b*NSB
           # indirect effect (a*b)
             ab := a*b
           # total effect
             total := c + (a*b)
         '

set.seed(123)
fit = sem(model, data = MED, se = "bootstrap", bootstrap = 1000)
summary(fit, fit.measures = TRUE, standardized = TRUE)
## lavaan 0.6-20 ended normally after 1 iteration
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                         5
## 
##   Number of observations                           191
## 
## Model Test User Model:
##                                                       
##   Test statistic                                 0.000
##   Degrees of freedom                                 0
## 
## Model Test Baseline Model:
## 
##   Test statistic                               176.345
##   Degrees of freedom                                 3
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    1.000
##   Tucker-Lewis Index (TLI)                       1.000
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -1353.544
##   Loglikelihood unrestricted model (H1)      -1353.544
##                                                       
##   Akaike (AIC)                                2717.088
##   Bayesian (BIC)                              2733.349
##   Sample-size adjusted Bayesian (SABIC)       2717.511
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.000
##   90 Percent confidence interval - lower         0.000
##   90 Percent confidence interval - upper         0.000
##   P-value H_0: RMSEA <= 0.050                       NA
##   P-value H_0: RMSEA >= 0.080                       NA
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.000
## 
## Parameter Estimates:
## 
##   Standard errors                            Bootstrap
##   Number of requested bootstrap draws             1000
##   Number of successful bootstrap draws             998
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   BDI ~                                                                 
##     ABK        (c)    0.114    0.088    1.297    0.195    0.114    0.083
##   NSB ~                                                                 
##     ABK        (a)    0.507    0.082    6.159    0.000    0.507    0.420
##   BDI ~                                                                 
##     NSB        (b)    0.771    0.078    9.866    0.000    0.771    0.681
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .BDI              60.722    6.349    9.564    0.000   60.722    0.482
##    .NSB              80.736    7.804   10.345    0.000   80.736    0.823
## 
## Defined Parameters:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##     ab                0.391    0.074    5.289    0.000    0.391    0.286
##     total             0.505    0.090    5.581    0.000    0.505    0.369

Der Zusammenhang zwischen Abhängigkeitskognitionen und der Ausprägung der Depressivität wird vollständig durch die negative Selbstbewertung vermittelt. Der direkte Effekte von Abhängigkeitskognitionen auf Depressivität ist nicht signifikant, während sowohl der Effekt von Abhängigkeitskognitionen auf die negative Selbstbewertung als auch der Effekt der negativen Selbstbewertung auf die Depressivität signifikant sind.